#source("/var/www/psych/dataR/runOnStart.R")
processData<-function() {
    #commented this out because it led to errors...look in apache err logs for errors
    #RApacheOutputErrors(status=TRUE,prefix='&lt!--\n',suffix='--&gt\n')
    
    #also helpful: sudo /etc/init.d/apache2 restart
    # and: tail -f  /var/log/apache2/error.log

    #todo: get dynamic directory
    dataRdir<-"/var/www/psych/dataR/"

    #d<-read.table("c:/temp/data216.185.22.90.tab")
    fname<-GET$filename    
    d<-read.table(file=paste(dataRdir,"data",fname,".tab",sep=""),sep="\t")
    info<-read.table(file=paste(dataRdir,"info",fname,".tab",sep=""),sep="\t",as.is=T)
    prog.goal<-F
    if(ncol(d)==3) { #mean value
        names(d)<-c("value","period","treated") 
    } else { #prog to goal
        prog.goal<-T
        names(d)<-c("value","n.size","g.days","g.value","period","treated") 
    }


    #dst for dataset, wt for wt, nss for n-sizes, wnss for weighted nss, mns for means, vrs for variances
    dst<-wt<-list(4)
    nss<-wnss<-mns<-vrs<-sds<-rep(NA,4)
    dst[[1]]<-matrix(d$value[d$treated==1 & d$period<0],ncol=1) #bef treat
    dst[[2]]<-matrix(d$value[d$treated==1 & d$period>=0],ncol=1) #aft treat
    dst[[3]]<-matrix(d$value[d$treated==0 & d$period<0],ncol=1) #bef control
    dst[[4]]<-matrix(d$value[d$treated==0 & d$period>=0],ncol=1) #aft control

    wt[[1]]<-exp(d$period[d$treated==1 & d$period<0]+1)
    wt[[2]]<-exp(-1*(d$period[d$treated==1 & d$period>=0]))
    wt[[3]]<-exp(d$period[d$treated==0 & d$period<0]+1)
    wt[[4]]<-exp(-1*(d$period[d$treated==0 & d$period>=0]))


    if (prog.goal) {
        control.goal.value<-d$g.value[d$treated==0][1]
        control.goal.days<-d$g.days[d$treated==0][1]
        treatment.goal.value<-d$g.value[d$treated==1][1]
        treatment.goal.days<-d$g.days[d$treated==1][1]
    
        #get the daily n-sizes
        daily.ns<-list(4)
        daily.ns[[1]]<-d$n.size[d$treated==1 & d$period<0] #bef treat
        daily.ns[[2]]<-d$n.size[d$treated==1 & d$period>=0] #aft treat
        daily.ns[[3]]<-d$n.size[d$treated==0 & d$period<0] #bef control
        daily.ns[[4]]<-d$n.size[d$treated==0 & d$period>=0] #aft control

        #fix both the DV and the weight
        for(i in 1:4) {
            dst[[i]]<-dst[[i]]/ifelse(i<3,(control.goal.value/control.goal.days),(treatment.goal.value/treatment.goal.days))
            wt[[i]]<-wt[[i]]*sqrt(daily.ns[[i]])
        }

    }

    for(i in 1:4) {
        if (!prog.goal) {
            nss[i]<-length(dst[[i]])
        } else {
            nss[i]<-sum(daily.ns[[i]])
        }
        wnss[i]<-sum(wt[[i]])
        mns[i]<-weighted.mean(dst[[i]],wt[[i]])
        if (nss[i]>1) {
            vrs[i]<-cov.wt(dst[[i]],wt[[i]])$cov[1,1]/mean(nss[i],wnss[i]) #not sure whether to use n size or weight size
            sds[i]<-sqrt(vrs[i])
        }
    }
    if(min(nss)>50 && min(wnss)>ifelse(prog.goal,5,20)) {
        treatment.effect.m<-(mns[2]-mns[1])-(mns[4]-mns[3])
        treatment.effect.v<-sum(vrs)
        treatment.effect.sd<-sqrt(treatment.effect.v)
        prob.t.e.not.zero<- pnorm(0,m=treatment.effect.m,sd=treatment.effect.sd)
        if (treatment.effect.m>0) {
            prob.t.e.not.zero<-(1 - prob.t.e.not.zero)
        }
        enough.data<-T
    } else {
        treatment.effect.m<-(mns[2]-mns[1])-(mns[4]-mns[3])
        prob.t.e.above.zero<-treatment.effect.sd<-treatment.effect.v<-NA
        enough.data<-F
    }
    
    #graph 1: overall treatment effect
    #to do: better color
    
    all.height<-380
    wider.width<-all.height
    mid.width<-wider.width*6/8
    thinner.width<-wider.width/2
    
    jpeg(filename=paste(dataRdir,"overall",fname,".jpg",sep=""),width=thinner.width,height=all.height)
    standard.par()
    if (enough.data) {
        ylim.most<-treatment.effect.m+treatment.effect.sd*3*c(-1,1) #make sure ylim covers sd
        ylim.good<-c(min(c(0,ylim.most[1])),max(c(0,ylim.most[2])))   #but also make sure it touches 0
        overall.color<-get.bar.col(treatment.effect.m,treatment.effect.sd)
        bp.x<-barplot(treatment.effect.m,ylim=ylim.good,col=overall.color)
        points(bp.x,treatment.effect.m,pch=19,cex=1.25)
        lines(rep(bp.x,2),treatment.effect.m+treatment.effect.sd*1*c(-1,1),lwd=2)
        lines(rep(bp.x,2),treatment.effect.m+treatment.effect.sd*1.96*c(-1,1),lwd=1,lty=2)
        abline(h=0)
    } else {
        create.blank.plot(F)
    }
    try(dev.off())

    #graph 2: bar plot of diff-in-diff
    jpeg(filename=paste(dataRdir,"detail",fname,".jpg",sep=""),width=mid.width,height=all.height)
    standard.par()
    treatment.region<-info[1,2]
    treatment.region<-paste(substring(treatment.region,1,10),ifelse(nchar(treatment.region)>10,".",""),sep="")
    if (enough.data) {
        bp.x<-barplot(height=matrix(mns,byrow=F,ncol=2),beside=T,names.arg=c(treatment.region,"Nearby"),cex.names=1.2,ylim=c(min(c(mns,0))*1.2,max(c(mns,0))*1.2))
        text(bp.x[1,],mns[c(1,3)],label="Before",pos=3)
        text(bp.x[2,],mns[c(2,4)],label="After",pos=3)
        se.col<-"gray"
        for(i in 1:4) {
            #points(bp.x[i],mns[i],pch=19,cex=1.25,col=se.col)
            lines(rep(bp.x[i],2),mns[i]+sds[i]*1*c(-1,1),lwd=2,col=se.col)
            lines(rep(bp.x[i],2),mns[i]+sds[i]*1.96*c(-1,1),lwd=1,lty=2,col=se.col)
        }
    } else {
        create.blank.plot(F)
    }
    try(dev.off())


    #graph 3: line plot over time.
    #first, figure out when we're out of data
    jpeg(filename=paste(dataRdir,"time",fname,".jpg",sep=""),width=wider.width,height=all.height)
    standard.par()
    if (enough.data) {
        max.period<-14
        min.daily.data<-30
        min.days<-5
        good.days <- 0
        daily.n.t<-daily.n.c<-rep(NA,max.period+1)
        for(i in 0:max.period) {
            daily.n.c[i+1]<-sum(d$treated==0 & d$period==i)
            daily.n.t[i+1]<-sum(d$treated==1 & d$period==i)
            #good.days = number of consecutive days with good data
            if ((good.days>=(i-1)) && (daily.n.c[i+1]>=min.daily.data) && (daily.n.t[i+1]>=min.daily.data)) {
                good.days<-i
            }
        }
        
        before.effect<-mns[4]-mns[3]
        enough.time.data<-(good.days>=min.days)
        if(enough.time.data) {
            x.s<-as.array(0:good.days)
            y.s<-apply(x.s,1,function(day) {get.day.mean.effect(d,before.effect,day)})
            daily.df<-data.frame(y=y.s,x=x.s)
            l.s.o<-loess(y~x,daily.df,span=1.1)
            l.s<-l.s.o$fitted
            #plot(x.s,l.s,type="l",xlab="Days",ylab="",ylim=c(min(c(0,1.2*l.s)),max(c(0,1.2*l.s))),col=overall.color)
            x.s.2<-as.array(seq(from=0,to=good.days,by=.1))
            l.s.2<-predict(l.s.o,newdata=x.s.2)
            plot(x.s.2,l.s.2,type="l",xlab="Days",ylab="",ylim=c(min(c(0,1.2*l.s)),max(c(0,1.2*l.s))),col=overall.color)
            abline(h=treatment.effect.m,lty=3,col=overall.color)
            text(x=min(x.s.2),y=treatment.effect.m,label=("Overall effect"),col=overall.color,adj=c(0,-.5))
            abline(h=0,lty=2)
        } else {
            create.blank.plot(T)
        }
    } else {
        create.blank.plot(T)
    }
    try(dev.off())
    alt.text.1<-paste("mean:",signif(treatment.effect.m,2),"Total n:",signif(sum(nss)),
        "Weighted treated n:",signif(wnss[4],1),sep=" ")
    
    alt.text.2<-paste("means bc ac bt at:",paste(signif(mns,2),collapse=" "),
        "st errs",paste(signif(sds,2),collapse=" "),
        "n sizes",paste(nss,collapse=" "),
        "and wtd n sizes",paste(signif(wnss,1),collapse=" "),
        sep=" ")

    if (enough.data) {
        if(enough.time.data) {
            alt.text.3<-paste("y.s",paste(apply(as.array(y.s),1,function(y) {signif(y,2)}),collapse=" "),
            "loess",paste(apply(as.array(l.s),1,function(y) {signif(y,2)}),collapse=" "),sep=" ")
        } else {
            alt.text.3<-paste("gdays",good.days,"n sizes t",paste(daily.n.t,collapse=" "),"n sizes c",paste(daily.n.c,collapse=" "),sep=" ")
        }
    } else {
        alt.text.3<-paste("NA")
    }

    setContentType("text/html")
    cat('<HTML><body>')
    t.e.str<-ifelse(prog.goal,to.pct.str(treatment.effect.m),round(treatment.effect.m,1))
    cat(paste(" Processed!",ifelse(enough.data,t.e.str,"XX"),
                ifelse(enough.data,round(prob.t.e.not.zero,2),"XX"),alt.text.1,alt.text.2,alt.text.3,sep="|"))
    cat('|</body>')
    cat("</BODY></HTML>\n")
    return(1)
}

#to do: make better
get.bar.col<-function(t.m,t.sd) {
    ifelse(t.m>0,"blue","red")
}

#could add standard error here as well
get.day.mean.effect<-function(d,before.effect,day) {
    after.effect<-mean(d$value[d$treated==1 & d$period==day]) - mean(d$value[d$treated==0 & d$period==day])
    diff.in.diff<-after.effect - before.effect
    return(diff.in.diff)
}

create.blank.plot<-function(wide.plot=F) {
    plot(0,0,type="n",axes=F,xlab="",ylab="")
    par(family="HersheySans",font=2)
    if(wide.plot) {
        text(0,0,label="Not enough data",cex=2)
    } else {
        text(0,.12*c(1,0,-1),label=c("Not","enough","data"),cex=2)
    }
        
}

standard.par<-function() {
    par(mar=c(3,2,2,2))
}

to.pct.str<-function(x) {
    return(paste(round(x*100,0),"%",sep=""))
}
